---
title: "Descriptive analysis of environmental protests"
author: "Endre Borbáth"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
    html_document:
        toc: true
        toc_depth: 4
---

```{r setup, include=FALSE, echo=FALSE}
rm(list = ls())

library(tidyr)
library(dplyr)
library(lubridate)
library(ggplot2)
library(forcats)
library(ggrepel)
library(cowplot)
library(countrycode)
library(corrplot)
library(lmtest)
library(sandwich)
library(texreg)
library(broom)
library(DescTools)
library(prais)
library(fmsb)
library(scales)

options(scipen = 999)

set.seed(1991)

path <- "/Users/eborbath/Library/CloudStorage/Dropbox-WZB/Endre Borbath/Papers/Environmental_protests/"

# load(file = paste0(path, "data/PEA_24042023.RData"))

PEA <- read.csv(paste0(path, "data/PEA_2000_2021/PEA_integrated_2000_2021_V2_cleaner_longcov.csv"))

PEA <- PEA %>%
  select(
        country_name, year, date, weighted_event, weighted_part_all, contains("issue_"), contains("actor_"),
        action_form
    ) %>%
    mutate(
        env_events = ifelse(issue_environment == "yes", weighted_event, 0),
        env_particip = ifelse(issue_environment == "yes", weighted_part_all, 0)
    ) %>%
    mutate(country_name = ifelse(country_name == "northern ireland", "united kingdom", country_name)) %>%
    mutate(country_iso = countryname(country_name, "iso2c")) %>%
  select(-covidissue_pro, -covidissue_anti, -actor_string1, -actor_string2) %>% 
  mutate(date2=ifelse(year==2021, date, as.character(NA))) %>% 
  mutate(date=ifelse(year==2021, as.character(NA), date)) %>% 
  mutate(date=ymd(date)) %>% 
  mutate(date2=dmy(date2)) %>% 
  mutate(date=ifelse(year==2021, date2, date)) %>% 
  mutate(date=as.Date(date)) %>% 
  select(-date2) %>%
  filter(!(country_iso %in% c("LU", "MT", "IS"))) %>% 
  mutate(month=month(date)) %>% 
  mutate(month = case_when(
    nchar(as.character(month)) == 1 ~ paste0("0", month),
    TRUE ~ as.character(month)
  )) %>%
  mutate(ym_character=paste(year, month, sep = "-")) %>% 
  mutate(ym = as.Date(paste0(ym_character, "-01"))) %>% 
  mutate(ym = factor(format(ym, "%Y-%m"), levels = unique(format(sort(ym), "%Y-%m"))))

# table(PEA$env_particip[PEA$ym_character=="2019-09" & PEA$env_particip>50000])
# View(PEA[PEA$ym_character=="2019-09" & PEA$env_particip>50000, ])

# sum(PEA$env_particip[PEA$ym_character=="2002-01"])
  
# View(PEA[PEA$ym_character=="2002-01",])
# table(PEA$issue_environment[PEA$country_iso=="IS"])


freq_table <- function(df, col_names) {
    lapply(col_names, function(x) {
        table(df[[x]])
    })
}

miss_table <- function(missing_var) {
    miss_df <- PEA_yearly %>%
        group_by(country_iso, year) %>%
        mutate(number_of_missings = sum(is.na(.data[[missing_var]]))) %>%
        ungroup() %>%
        mutate(number_of_missings = ifelse(number_of_missings == 1, "Miss.", as.character(NA))) %>%
        filter(number_of_missings == "Miss.") %>%
        dplyr::select(country_iso, year, number_of_missings) %>%
        pivot_wider(
            names_from = country_iso,
            values_from = number_of_missings
        )

    print(miss_df, n = Inf)
}

class_table <- function(df, col_names) {
    lapply(col_names, function(x) {
        class(df[[x]])
    })
}

range01 <- function(x) {
    (x - min(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE))
}

process_model <- function(model) {
  model_name <- deparse(substitute(model))
  mod1_res <- tidy(model)
  mod1_ci <- as.data.frame(confint(model))
  mod1_ci$term <- rownames(mod1_ci)
  mod1_res <- merge(mod1_ci, mod1_res, by = "term", all = TRUE)
  mod1_ci <- as.data.frame(confint(model, level = 0.84))
  mod1_ci$term <- rownames(mod1_ci)
  mod1_res <- merge(mod1_ci, mod1_res, by = "term", all = TRUE)
  mod1_res$model_name <- model_name
  return(mod1_res)
}

```

## Overall in the dataset

```{r, echo=FALSE, message=FALSE, fig.width=8, fig.height=6}
# using ggplot2 create a scatter plot of the total of weighted_event by year

plot_dat <- PEA %>%
    select(env_events, env_particip) %>%
    summarize_all(~ sum(., na.rm = TRUE)) 

plot_dat <- PEA %>%
    # filter(action_form!="petitions, symbolic actions") %>% 
    select(ym_character, env_events, env_particip, weighted_event) %>%
    group_by(ym_character) %>%
    mutate_at(vars(starts_with("env_"), weighted_event), ~ sum(., na.rm = TRUE)) %>%
    ungroup() %>%
    unique() %>%
    # mutate(env_particip = env_particip / env_events) %>%
    mutate(env_event_share = env_events / weighted_event) %>%
    select(-weighted_event) %>%
    pivot_longer(cols = c(env_events, env_particip, env_event_share), names_to = "variable", values_to = "value") %>%
    mutate(variable = fct_recode(variable, "Number of events" = "env_events", "Number of participants" = "env_particip", "Share of events" = "env_event_share")) %>%
    mutate(variable = factor(variable, levels = c("Number of events", "Number of participants", "Share of events")))

# Convert year_month to Date object
plot_dat$ym <- as.Date(paste0(plot_dat$ym_character, "-01"))

# to avoid the two outliers of 2002 January (large petitions in Austria against the Nuclear power plant in Temelin, FPO sponsored.) and 2019 September (climate strikes in Germany and Italy), I will remove them

plot_dat <- plot_dat %>%
    mutate(value = ifelse((ym_character == "2002-01" | ym_character == "2019-09") & variable=="Number of participants", NA, value))

ggplot(plot_dat, aes(x = ym, y = value)) +
    facet_wrap(~variable, scales = "free_y", nrow = 2) +
    geom_point(alpha=0.2) +
    geom_smooth(method = "loess", se = FALSE, color = "black", size=0.75, span=0.75) +
      scale_x_date(breaks = seq(min(plot_dat$ym), max(plot_dat$ym), by = "3 year"), labels = date_format("%Y")) +
  theme_linedraw() +
    xlab("") +
    ylab("")

ggsave(
    path = paste0(path, "draft/July_2024/Figures/"),
    filename = "PEA_trends_overall.png",
    # scale = ,
    width = 6, height = 6,
    units = "in", dpi = 600
)
```


## By country - facets {.tabset}

### Events

```{r, echo=FALSE, message=FALSE, fig.width=8, fig.height=6}
# using ggplot2 create a scatter plot of the total of weighted_event by country

plot_dat <- PEA %>%
    select(country_iso, ym_character, env_events) %>%
    group_by(country_iso, ym_character) %>%
    mutate_at(vars(starts_with("env_")), ~ sum(., na.rm = TRUE)) %>%
    ungroup() %>%
    unique() %>%
    group_by(country_iso) %>%
    mutate(indic_mean = mean(env_events, na.rm = TRUE)) %>%
    ungroup() %>%
    mutate(country_iso = countrycode(country_iso, "iso2c", "country.name")) %>%
    mutate(country_iso = fct_reorder(country_iso, indic_mean, .desc = TRUE)) %>%
    arrange(indic_mean, ym_character)

means <- plot_dat %>%
    select(country_iso, indic_mean) %>%
    unique()

plot_dat$ym <- as.Date(paste0(plot_dat$ym_character, "-01"))

ggplot(plot_dat, aes(x = ym, y = env_events)) +
    geom_point(alpha=0.1) +
    geom_smooth(method = "loess", se = FALSE, color = "black", size=0.75, span=0.25) +
    scale_x_date(breaks = seq(min(plot_dat$ym), max(plot_dat$ym), by = "7 year"), labels = date_format("%Y")) +
    theme_linedraw() +
    geom_vline(aes(xintercept = as.Date("2019-01-01")), color = "black", linetype = "dashed") +
    geom_hline(
        data = means, aes(yintercept = indic_mean),
        color = "black", linetype = "dashed"
    ) +
    facet_wrap(~country_iso, scales = "free_y", nrow = 4) +
    scale_y_continuous(
        labels = scales::number_format(accuracy = 0.1)
    ) +
    labs(x = "Year", y = "Total number of events per year") +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.title.x = element_blank()
    )

ggsave(
    path = paste0(path, "draft/July_2024/Appendix/Figures/"),
    filename = "PEA_events_by_country_facets.png",
    scale = 1.4,
    width = 8, height = 6,
    units = "in", dpi = 300
)
```

### Participants

```{r, echo=FALSE, message=FALSE, fig.width=8, fig.height=6}
# using ggplot2 create a scatter plot of the total of weighted_event by country

plot_dat <- PEA %>%
    select(country_iso, ym_character, env_particip, env_events) %>%
    group_by(country_iso, ym_character) %>%
    mutate_at(vars(starts_with("env_")), ~ sum(., na.rm = TRUE)) %>%
    ungroup() %>%
    unique() %>%
    group_by(country_iso) %>%
    mutate(indic_mean = mean(env_events, na.rm = TRUE)) %>%
    ungroup() %>%
    select(-env_events) %>%
    mutate(country_iso = countrycode(country_iso, "iso2c", "country.name")) %>%
    mutate(country_iso = fct_reorder(country_iso, indic_mean, .desc = TRUE)) %>%
    arrange(indic_mean, ym_character)

means <- plot_dat %>%
    select(country_iso, env_particip) %>%
    group_by(country_iso) %>% 
    mutate(indic_mean = mean(env_particip, na.rm = TRUE)) %>%
    ungroup() %>%
    unique()

plot_dat$ym <- as.Date(paste0(plot_dat$ym_character, "-01"))

ggplot(plot_dat, aes(x = ym, y = env_particip/1000)) +
    geom_point(alpha=0.1) +
    geom_smooth(method = "loess", se = FALSE, color = "black", size=0.75, span=0.25) +
    scale_x_date(breaks = seq(min(plot_dat$ym), max(plot_dat$ym), by = "7 year"), labels = date_format("%Y")) +
    theme_linedraw() +
    geom_vline(aes(xintercept = as.Date("2019-01-01")), color = "black", linetype = "dashed") +
    geom_hline(
        data = means, aes(yintercept = indic_mean/1000),
        color = "black", linetype = "dashed"
    ) +
    facet_wrap(~country_iso, scales = "free_y", nrow = 4) +
    scale_y_continuous(
        labels = scales::number_format(accuracy = 0.1)
    ) +
    labs(x = "Year", y = "Total Number of participants in env. events per year (in thousands)") +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.title.x = element_blank()
    )

ggsave(
    path = paste0(path, "draft/July_2024/Appendix/Figures/"),
    filename = "PEA_participants_by_country_facets.png",
    scale = 1.4,
    width = 8, height = 6,
    units = "in", dpi = 300
)
```

### Share of environmental events

```{r, echo=FALSE, message=FALSE, fig.width=8, fig.height=6}
# using ggplot2 create a scatter plot of the total of weighted_event by country

plot_dat <- PEA %>%
    select(country_iso, ym_character, env_events, weighted_event) %>%
    group_by(country_iso, ym_character) %>%
    mutate(share_event = sum(env_events, na.rm = TRUE) /
        sum(weighted_event, na.rm = TRUE)) %>%
    mutate(env_events = sum(env_events, na.rm = TRUE)) %>%
    ungroup() %>%
    select(country_iso, ym_character, share_event, env_events) %>%
    unique() %>%
    group_by(country_iso) %>%
    mutate(indic_mean = mean(env_events, na.rm = TRUE)) %>%
    ungroup() %>%
    mutate(country_iso = countrycode(country_iso, "iso2c", "country.name")) %>%
    mutate(country_iso = fct_reorder(country_iso, indic_mean, .desc = TRUE)) %>%
    arrange(indic_mean, ym_character) %>%
    mutate(indic_mean = mean(share_event, na.rm = TRUE), .by = "country_iso")

means <- plot_dat %>%
    select(country_iso, share_event) %>%
    group_by(country_iso) %>% 
    mutate(indic_mean = mean(share_event, na.rm = TRUE)) %>%
    ungroup() %>%
    unique()

plot_dat$ym <- as.Date(paste0(plot_dat$ym_character, "-01"))

ggplot(plot_dat, aes(x = ym, y = share_event)) +
    geom_point(alpha=0.1) +
    geom_smooth(method = "loess", se = FALSE, color = "black", size=0.75, span=0.25) +
    scale_x_date(breaks = seq(min(plot_dat$ym), max(plot_dat$ym), by = "7 year"), labels = date_format("%Y")) +
    theme_linedraw() +
    geom_vline(aes(xintercept = as.Date("2019-01-01")), color = "black", linetype = "dashed") +
    geom_hline(
        data = means, aes(yintercept = indic_mean),
        color = "black", linetype = "dashed"
    ) +
    facet_wrap(~country_iso, scales = "free_y", nrow = 4) +
    scale_y_continuous(
        labels = scales::number_format(accuracy = 0.1)
    ) +
    labs(x = "Year", y = "Share of environmental events from all protest events") +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.title.x = element_blank()
    )

ggsave(
    path = paste0(path, "draft/July_2024/Appendix/Figures/"),
    filename = "PEA_share_events_by_country.png",
    scale = 1.4,
    width = 8, height = 6,
    units = "in", dpi = 300
)
```


# Cross-sectional analysis of environmental events

```{r, echo=FALSE, message=FALSE, fig.width=8, fig.height=6}
# using ggplot2 create a scatter plot of the total of weighted_event by year

plot_dat <- PEA %>%
    select(country_iso, env_events, env_particip, weighted_event) %>%
    group_by(country_iso) %>%
    mutate_at(vars(starts_with("env_")), ~ sum(., na.rm = TRUE)) %>%
    mutate(weighted_event = sum(weighted_event, na.rm = TRUE)) %>%
    ungroup() %>%
    unique() %>%
    mutate(env_event_share = env_events / weighted_event) %>%
    # mutate(new_env_particip = env_particip / env_events) %>%
    # mutate(new_env_particip = ifelse(env_particip == 0 & env_events == 0, 0, new_env_particip)) %>%
    # mutate(env_particip = new_env_particip) %>%
    # select(-new_env_particip) %>%
    ## mutate(env_particip = env_particip / 100000) %>%
    mutate(country_iso = fct_reorder(country_iso, env_events, .desc = TRUE)) %>%
    mutate(dropped_countries = ifelse(env_events <= 20, "Dropped countries", "Kept countries")) %>%
    mutate(dropped_countries = factor(dropped_countries, levels = c("Kept countries", "Dropped countries"))) %>%
    arrange(env_events) %>%
    select(-weighted_event) %>%
    pivot_longer(cols = c(env_events, env_particip, env_event_share), names_to = "variable", values_to = "value") %>%
    mutate(variable = fct_recode(variable, "Number of events" = "env_events", "Number of participants" = "env_particip", "Share of events" = "env_event_share")) %>%
    mutate(variable = factor(variable, levels = c("Number of events", "Number of participants", "Share of events")))

ggplot(plot_dat, aes(x = country_iso, y = value, fill = dropped_countries)) + #
    geom_bar(stat = "identity", color = "black") +
    facet_wrap(~variable, scales = "free_y", nrow = 1) +
    theme_bw() +
    scale_fill_manual(values = c("gray40", "gray90")) +
    xlab("") +
    ylab("") +
    theme(
        legend.position = "bottom", legend.title = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
    )

ggsave(
    path = paste0(path, "draft/July_2024/Appendix/Figures/"),
    filename = "PEA_country_means.png",
    scale = 1.6,
    width = 6, height = 4,
    units = "in", dpi = 300
)
```

```{r, echo=FALSE, message=FALSE, fig.width=8, fig.height=6}
load(file = paste0(path, "data/PEA_yearly_16052024.RData"))

cor_dat <- PEA_yearly %>%
    mutate(green_in_parl = ifelse(green_seat_share > 0,
        1, 0
    )) %>%
    mutate_at(vars(green_in_gov, green_in_parl), ~ as.numeric(.)) %>%
    select(
        total_env_protest, total_env_part, share_env_protest,
        econ_indic, EPI, climate_change_index,
        env_health_index, eco_vitality_index, nr_disasters,
        gov_env_cmp, env_prot_sys,
        green_in_gov, green_in_parl,
        share_organized,
        NGO_IUCN, share_demo
    ) %>%
    mutate_at(vars(everything()), ~ as.numeric(.))

M <- cor(cor_dat, use = "pairwise.complete.obs")
head(round(M, 2))

rownames(M) <- c(
    "Number of env. events",
    "Part. in env. events",
    "Share of env. events",
    "Economic performance",
    "Environmental performance",
    "Climate change",
    "Environmental health",
    "Ecosystem vitality",
    "Number of natural disasters",
    "Gov. position on env. protection",
    "Parl. position on env. protection",
    "Green party in gov.",
    "Green party in parl.",
    "Share of organized protests",
    "Number of env. NGOs",
    "Share of demonstrations"
)

colnames(M) <- c(
    "Number of env. events",
    "Part. in env. events",
    "Share of env. events",
    "Economic performance",
    "Environmental performance",
    "Climate change",
    "Environmental health",
    "Ecosystem vitality",
    "Number of natural disasters",
    "Gov. position on env. protection",
    "Parl. position on env. protection",
    "Green party in gov.",
    "Green party in parl.",
    "Share of organized protests",
    "Number of env. NGOs",
    "Share of demonstrations"
)

col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))

png(
    filename = paste0(path, "draft/July_2024/Appendix/Figures/corplot_yearly.png"),
    width = 7, height = 7, res = 500, units = "cm", pointsize = 3
)
corrplot(M,
    method = "color", col = col(200),
    type = "upper",
    addCoef.col = "black", # Add coefficient of correlation
    tl.col = "black", tl.srt = 45, # Text label color and rotation
    diag = FALSE,
    mar = c(rep(0, 4)),
    number.cex = 0.8
)
dev.off()
```

```{r, echo=FALSE, message=FALSE, fig.width=8, fig.height=6}

# dput(colnames(PEA)[grepl("actor", colnames(PEA))])

PEA_reg <- PEA %>%
    mutate(only_env = ifelse(
        issue_environment == "yes" & issue_econ_private == "no" & issue_econ_public == "no" &
            issue_cult_lib == "no" & issue_political == "no" & issue_regional == "no" & issue_cult_cons == "no" &
            issue_xeno == "no" & issue_other == "no" & issue_anti_europe == "no" & issue_anti_immig == "no" &
            issue_pro_europe == "no" & issue_pro_immig == "no" & issue_corona_anti_restric == "no" & issue_corona_pro_restric == "no" &
            issue_education == "no" & issue_healthcare == "no",
        as.numeric(1), as.numeric(0)
    )) %>%
    mutate(
        issue_econ = ifelse(issue_econ_private == "yes" | issue_econ_public == "yes", "yes", "no"),
        issue_cult = ifelse(issue_cult_lib == "yes" | issue_cult_cons == "yes" | issue_xeno == "yes" |
                            issue_anti_immig == "yes" | issue_pro_immig == "yes" |
                            issue_corona_anti_restric == "yes" | issue_corona_pro_restric == "yes", "yes", "no"),
        issue_political = ifelse(issue_political == "yes" | issue_regional == "yes" |
                                 issue_anti_europe == "yes" | issue_pro_europe == "yes", "yes", "no"),
        issue_other = ifelse(issue_other == "yes" | issue_dontknow  == "yes" | issue_missing == "yes" | 
                               issue_education == "yes" | issue_healthcare == "yes", "yes", "no")
    ) %>%
    mutate_at(
        vars(
            issue_econ, issue_cult, issue_political, issue_other, issue_environment),
        ~ as.numeric(case_when(
            . == "yes" ~ 1,
            . == "no" ~ 0,
            TRUE ~ 0
        ))
    ) %>%
    mutate_at(vars(issue_econ, issue_cult, issue_political, issue_other), ~ as.factor(.)) %>%
    mutate(form = case_when(
        action_form %in% c("confrontations, blockades") ~ "Confrontative protest",
        action_form %in% c("demonstrations") ~ "Demonstrations",
        action_form %in% c("petitions, symbolic actions") ~ "Petitions & symbolic actions",
        action_form %in% c("violent protest") ~ "Violent protest",
        action_form %in% c("strikes") ~ "Strikes",
        action_form %in% c("None", "other protest") ~ "Other forms"
    )) %>%
    mutate(form = factor(form, levels = c("Demonstrations", 
                                          "Petitions & symbolic actions",
                                          "Strikes",
                                          "Confrontative protest", 
                                          "Violent protest", 
                                          "Other forms"))) %>%
    mutate(
        demonstrations = ifelse(form == "Demonstrations", 1, 0),
        petitions = ifelse(form=="Petitions & symbolic actions", 1, 0),
        strikes = ifelse(form=="Strikes", 1, 0),
        confrontative = ifelse(form == "Confrontative protest", 1, 0), 
        violent = ifelse(form == "Violent protest", 1, 0), 
        other = ifelse(form == "Other forms", 1, 0)
    ) %>%
    mutate_at(vars(demonstrations, petitions, strikes, 
                   confrontative, violent, other), ~ as.factor(.)) %>%
    mutate(no_org_actor = ifelse((actor_party_left == "no" & actor_party_right == "no" & actor_party_unknown == "no" &
        actor_union_private == "no" & actor_union_public == "no" & actor_union_both == "no" &
        actor_union_unknown == "no" &
        actor_party_far_left == "no" & actor_party_far_right == "no" & actor_party_main_left == "no" &
        actor_party_main_right == "no" & actor_prof_org == "no" & actor_non_prof_org == "no") &
        (actor_dontknow == "yes" | actor_missing == "yes" | actor_general_citizens == "yes" |
            actor_group_pens == "yes" | actor_group_stud == "yes" |
            actor_group_occup == "yes" | actor_group_other == "yes" | actor_other == "yes" |
            actor_group_lgbt == "yes" | actor_group_refugee == "yes" | actor_group_women == "yes"),
    as.numeric(1), as.numeric(0)
    )) %>%
    mutate(prof_org = ifelse(actor_prof_org == "yes", 1, 0)) %>%
    mutate(non_prof_org = ifelse(actor_non_prof_org == "yes", 1, 0)) %>%
    mutate(unions = ifelse(actor_union_private == "yes" | actor_union_public == "yes" | actor_union_both == "yes" |
        actor_union_unknown == "yes", 1, 0)) %>%
    mutate(parties = ifelse(actor_party_left == "yes" | actor_party_right == "yes" | actor_party_unknown == "yes" |
        actor_party_far_left == "yes" | actor_party_far_right == "yes" | actor_party_main_left == "yes" |
        actor_party_main_right == "yes", 1, 0)) %>%
    mutate_at(vars(prof_org, non_prof_org, unions, parties), ~ ifelse(is.na(.), 0, .)) %>%
    group_by(country_iso, year) %>%
    mutate(big_event = ifelse(weighted_part_all >= quantile(weighted_part_all, 0.80), 1, 0)) %>%
    mutate(big_event = as.factor(big_event)) %>% 
    mutate(year=as.factor(year)) %>% 
    mutate(weighted_part_all=range01(weighted_part_all)) 
# table(PEA_reg$issue_10, useNA = "always") # there is some missing, but all of them are accounted for by issue_dontknow=="yes" & issue_missing=="yes"

M1_fe <- glm(
    issue_environment ~ issue_econ + issue_cult + issue_political + 
      issue_other + petitions + strikes + confrontative + violent + other + 
    prof_org + non_prof_org + unions + parties + weighted_part_all + country_iso + year,
    data = PEA_reg, family = binomial(link = "logit")
)

M1_rob <- coeftest(M1_fe, vcovCL(M1_fe, cluster = list(PEA_reg$year, PEA_reg$country_iso)))

# dim(subset(PEA_reg, year %in% c(2019:2021)))
# dim(subset(PEA_reg, year %in% c(2000:2018)))


M2_fe <- glm(
    issue_environment ~ issue_econ + issue_cult + issue_political + 
      issue_other + petitions + strikes + confrontative + violent + other + 
    prof_org + non_prof_org + unions + parties + weighted_part_all + country_iso + year,
    data = subset(PEA_reg, year %in% c(2000:2018)), family = binomial(link = "logit")
)

M2_rob <- coeftest(M2_fe, vcovCL(M2_fe, cluster = list(subset(PEA_reg, year %in% c(2000:2018))$year, subset(PEA_reg, year %in% c(2000:2018))$country_iso)))

M3_fe <- glm(
    issue_environment ~ issue_econ + issue_cult + issue_political + 
      issue_other + petitions + strikes + confrontative + violent + other + 
    prof_org + non_prof_org + unions + parties + weighted_part_all + country_iso + year,
    data = subset(PEA_reg, year %in% c(2019:2021)), family = binomial(link = "logit")
)

M3_rob <- coeftest(M3_fe, vcovCL(M3_fe, cluster = list(subset(PEA_reg, year %in% c(2019:2021))$year, subset(PEA_reg, year %in% c(2019:2021))$country_iso)))

```

```{r, echo=FALSE, message=FALSE, fig.width=8, fig.height=6}
texreg(
    l = list(M1_rob, M2_rob, M3_rob), #doctype = FALSE, center = TRUE, siunitx = TRUE, use.packages = FALSE, booktabs = TRUE, float.pos = "h!", single.row = TRUE, no.margin = TRUE,
    caption = "Differences between environmental and other events",
   custom.model.names = c("Model 1" = "Full sample", "Model 2" = "2000-2018", "Model 3" = "2019-2021"),
   float.pos = "H",
    custom.gof.rows = list(
        "Country fixed effects" = rep("Yes", 3),
        "Year fixed effects" = rep("Yes", 3),
        "Adjusted R sq" = c(NagelkerkeR2(M1_fe)$R2, NagelkerkeR2(M2_fe)$R2, NagelkerkeR2(M3_fe)$R2),
        "Number of observations" = c(nobs(M1_fe), nobs(M2_fe), nobs(M3_fe))),
    omit.coef = "\\b(Intercept|year|country\\(\\w+\\))\\b",
    custom.coef.map = list(
        "issue_econ1" = "Economic",
        "issue_cult1" = "Cultural",
        "issue_political1" = "Political",
        "issue_other1" = "Other issue",
        "petitions1" = "Petitions \\& symbolic actions",
        "strikes1" = "Strikes",
        "confrontative1" = "Confrontative protest", 
        "violent1" = "Violent protest", 
        "other1" = "Other forms",
        "parties" = "Parties",
        "unions" = "Unions",
        "prof_org" = "Prof. organizations",
        "non_prof_org" = "Non-prof. organizations",
        "weighted_part_all" = "Participation rate"
    ),
    groups = list(
        "Combined issues (ref: env. alone)" = 1:4,
        "Form (ref: demonstrative)" = 5:7,
        "Organizers (ref: no organization)" = 8:11,
        "Size" = 12),
    include.deviance = FALSE,
    include.loglik = FALSE,
    file = paste0(path, "draft/July_2024/Appendix/Tables/events_clusteredSE.tex")
)
```



```{r, echo=FALSE, message=FALSE, fig.width=8, fig.height=6}

plot_dat <- process_model(M1_rob) %>%
    filter(!grepl("country_", term)) %>%
    filter(!grepl("Intercept", term)) %>%
    filter(!grepl("year", term)) %>%
    mutate(
        odds = exp(estimate),
        odds.outter.upper = exp(`2.5 %`),
        odds.outter.lower = exp(`97.5 %`),
        odds.inner.upper = exp(`8 %`),
        odds.inner.lower = exp(`92 %`)
    ) %>%
  # mutate_at(vars(contains("odds")), ~ifelse(.<1, -1/., .)) %>%
    mutate(term_orig = term) %>%
    mutate(term = case_when(
        term=="issue_econ1" ~ "Economic",
term=="issue_cult1" ~ "Cultural",
term=="issue_political1" ~ "Political",
term=="issue_other1" ~ "Other issue",
term=="petitions1" ~ "Petitions & symbolic actions",
term=="strikes1" ~ "Strikes",
term=="confrontative1" ~ "Confrontative protest", 
term=="violent1" ~ "Violent protest", 
term=="other1" ~ "Other forms",
term=="parties" ~ "Parties",
term=="unions" ~ "Unions",
term=="prof_org" ~ "Prof. organizations",
term=="non_prof_org" ~ "Non-prof. organizations",
term=="weighted_part_all" ~ "Participation rate",
        TRUE ~ term
    ))

titles <- plot_dat %>%
    filter(term_orig %in% c("issue_cult1", "violent1", "unions", "prof_org")) %>% # I just need some unique values
    mutate_at(vars(c(2:14)), ~NA)

titles$term <- c(
    "Combined issues (ref: env. alone)", "Form (ref: demonstrations)", "Organizers (ref: no organization)", "Size"
)


plot_dat <- bind_rows(plot_dat, titles)

plot_dat <- plot_dat %>%
    mutate(term = factor(term, levels = c(
         "Combined issues (ref: env. alone)", "Economic", "Cultural", "Political", "Other issue",
  "Form (ref: demonstrations)", "Petitions & symbolic actions", "Strikes", "Confrontative protest",  "Violent protest",  "Other forms",
  "Organizers (ref: no organization)", "Parties", "Unions", "Prof. organizations", "Non-prof. organizations",
  "Size", "Participation rate"
    ))) %>% 
  mutate(term=fct_rev(term))

ggplot(plot_dat, aes(y = term)) +
    geom_vline(xintercept = 1, linetype = "dashed", color = "gray60") +
    geom_point(aes(x = odds), size = 2) +
    geom_linerange(aes(
        xmin = odds.outter.lower, xmax = odds.outter.upper,
    ), size = 0.5, key_glyph = "path") +
    geom_linerange(aes(
        xmin = odds.inner.lower, xmax = odds.inner.upper,
    ), size = 1.1, key_glyph = "path") +
    xlab("Odds Ratio") +
    scale_color_brewer(name = "", palette = "Set2", guide = guide_legend(reverse = TRUE, nrow = 2, byrow = TRUE)) +
    scale_shape_discrete(name = "", guide = guide_legend(reverse = TRUE, nrow = 2, byrow = TRUE)) +
    theme_linedraw() +
    theme(
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.key.width = unit(1, "cm"),
        axis.title.y = element_blank(),
        axis.text.y = element_text(face = c(
             rep("plain", 1), "bold",
            rep("plain", 4), "bold",
            rep("plain", 5), "bold",
            rep("plain", 4), "bold"
        )),
        # axis.text.y = element_text(face = c(rep(2, 4))),
        axis.ticks.y = element_blank()
    )

ggsave(
    plot = last_plot(),
    filename = "coefplot_env_events.eps",
    path = paste0(path, "draft/July_2024/Figures/"),
    width = 7,
    height = 4,
    scale = 0.9,
    dpi = 400
)
```

```{r, echo=FALSE, message=FALSE, fig.width=8, fig.height=6}

plot_dat <- bind_rows(process_model(M2_rob), process_model(M3_rob))

plot_dat <- plot_dat %>%
    filter(!grepl("country_", term)) %>%
    filter(!grepl("Intercept", term)) %>%
    filter(!grepl("year", term)) %>%
    mutate(type=case_when(model_name=="M2_rob" ~ "2000-2018",
                          model_name=="M3_rob" ~ "2019-2021",
                          TRUE ~ as.character(NA))) %>% 
    mutate(
        odds = exp(estimate),
        odds.outter.upper = exp(`2.5 %`),
        odds.outter.lower = exp(`97.5 %`),
        odds.inner.upper = exp(`8 %`),
        odds.inner.lower = exp(`92 %`)
    ) %>%
    mutate(term_orig = term) %>%
    mutate(term = case_when(
        term=="issue_econ1" ~ "Economic",
term=="issue_cult1" ~ "Cultural",
term=="issue_political1" ~ "Political",
term=="issue_other1" ~ "Other issue",
term=="petitions1" ~ "Petitions & symbolic actions",
term=="strikes1" ~ "Strikes",
term=="confrontative1" ~ "Confrontative protest", 
term=="violent1" ~ "Violent protest", 
term=="other1" ~ "Other forms",
term=="parties" ~ "Parties",
term=="unions" ~ "Unions",
term=="prof_org" ~ "Prof. organizations",
term=="non_prof_org" ~ "Non-prof. organizations",
term=="weighted_part_all" ~ "Participation rate",
        TRUE ~ term
    ))

titles <- plot_dat %>%
    filter(term_orig %in% c("issue_cult1", "violent1", "unions", "prof_org")) %>% # I just need some unique values
    mutate_at(vars(c(2:14)), ~NA)

titles$term <- c(
    "Combined issues (ref: env. alone)", "Form (ref: demonstrations)", "Organizers  (ref: no organization)", "Size"
)


plot_dat <- bind_rows(plot_dat, titles)

plot_dat <- plot_dat %>%
    mutate(term = factor(term, levels = c(
         "Combined issues (ref: env. alone)", "Economic", "Cultural", "Political", "Other issue",
  "Form (ref: demonstrations)", "Petitions & symbolic actions", "Strikes", "Confrontative protest",  "Violent protest",  "Other forms",
  "Organizers  (ref: no organization)", "Parties", "Unions", "Prof. organizations", "Non-prof. organizations",
  "Size", "Participation rate"
    ))) %>% 
  mutate(term=fct_rev(term))

space_between_bars <- 0.65  
  
ggplot(plot_dat, aes(y = term, group=type)) +
    geom_vline(xintercept = 1, linetype = "dashed", color = "gray60") +
    geom_point(aes(x = odds, shape = type, color = type), position = position_dodge(width = space_between_bars), size = 2) +
    geom_linerange(aes(
        xmin = odds.outter.lower, xmax = odds.outter.upper, group = type, color = type
    ), size = 0.5, position = position_dodge(width = space_between_bars), key_glyph = "path") +
    geom_linerange(aes(
        xmin = odds.inner.lower, xmax = odds.inner.upper,  group = type, color = type
    ), size = 1.1, position = position_dodge(width = space_between_bars), key_glyph = "path") +
    xlab("Odds Ratio") +
    scale_color_brewer(name = "", palette = "Set2", na.translate = F,
                       guide = guide_legend(reverse = TRUE, nrow = 1, byrow = TRUE)) +
    scale_shape_discrete(name = "", na.translate = F,
                         guide = guide_legend(reverse = TRUE, nrow = 1, byrow = TRUE)) +
    theme_linedraw() +
    theme(
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.key.width = unit(1, "cm"),
        axis.title.y = element_blank(),
        axis.text.y = element_text(face = c(
             rep("plain", 1), "bold",
            rep("plain", 4), "bold",
            rep("plain", 5), "bold",
            rep("plain", 4), "bold"
        )),
        # axis.text.y = element_text(face = c(rep(2, 4))),
        axis.ticks.y = element_blank()
    )

ggsave(
    plot = last_plot(),
    filename = "coefplot_env_events_by_2019.png",
    path = paste0(path, "draft/July_2024/Appendix/Figures/"),
    width = 7,
    height = 6,
    scale = 0.9,
    dpi = 400
)
```

# Regression models - over time

```{r, echo=FALSE, message=FALSE, fig.width=8, fig.height=6}
subset_reg_data <- function(threshold) {
    reg_data <- PEA_yearly %>%
        mutate(green_in_parl = ifelse(green_seat_share > 0,
            1, 0
        )) %>%
        mutate_at(vars(green_in_gov, green_in_parl), ~ as.numeric(.)) %>%
        mutate(NGO_IUCN = ifelse(is.na(NGO_IUCN), 0, NGO_IUCN)) %>%
        ungroup() %>%
        mutate(country_iso = as.factor(country_iso)) %>%
        group_by(country_iso) %>%
        mutate(total_env_country = sum(total_env_protest, na.rm = TRUE)) %>%
        ungroup() %>%
        filter(total_env_country >= threshold)

    return(reg_data)
}
```




```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}

reg_data <- subset_reg_data(threshold = 20)

table(is.na(reg_data$year))

dput(sort(countrycode::countrycode(unique(PEA_yearly$country_iso), origin = "iso2c", destination = "country.name")))

reg_data <- reg_data %>%
    select(
        total_env_protest, total_env_part, share_env_protest,
        econ_indic,  EPI, climate_change_index,
        env_health_index, eco_vitality_index, nr_disasters, share_demo,
        gov_env_cmp, env_prot_sys, green_in_gov, green_in_parl,
        share_organized, NGO_IUCN,
        country_iso, year
    ) 

Amelia::missmap(reg_data)

reg_data <- reg_data %>%
    na.omit() %>%
    mutate_at(vars(
        total_env_protest, total_env_part, share_env_protest, econ_indic, 
        EPI, climate_change_index,
        env_health_index, eco_vitality_index, nr_disasters, share_demo,
        gov_env_cmp, env_prot_sys, green_in_gov, green_in_parl,
        share_organized, NGO_IUCN
    ), ~ scale(.))

rsqs <- as.data.frame(expand.grid(
    type = c("Grievances", "Opportunities", "Resources"),
    model = c("Events", "Participants", "Share events"),
    adjr2 = as.numeric(NA)
))
```



## Grievances

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
form_events <- as.formula(total_env_protest ~
    econ_indic +
    EPI +
    nr_disasters +
    share_demo +
    country_iso)

form_participants <- as.formula(total_env_part ~
    econ_indic +
    EPI +
    nr_disasters +
    share_demo +
    country_iso)

form_share_events <- as.formula(share_env_protest ~
    econ_indic +
    EPI +
    nr_disasters +
    share_demo +
    country_iso)


M1 <- prais_winsten(
    formula = form_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)

# M1_se <- sqrt(diag(vcovPC(M1)))
# coefci(M1, level = 0.95, vcov. = vcovPC(M1))

M1_rob <- coeftest(M1, vcov = vcovPC(M1))

M2 <- prais_winsten(
    formula = form_participants,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)

M2_rob <- coeftest(M2, vcov = vcovPC(M2))

M3 <- prais_winsten(
    formula = form_share_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)

M3_rob <- coeftest(M3, vcov = vcovPC(M3))

rsqs$adjr2[rsqs$type == "Grievances" & rsqs$model == "Events"] <- summary(M1)$adj.r.squared
rsqs$adjr2[rsqs$type == "Grievances" & rsqs$model == "Participants"] <- summary(M2)$adj.r.squared
rsqs$adjr2[rsqs$type == "Grievances" & rsqs$model == "Share events"] <- summary(M3)$adj.r.squared

# Print regression table
texreg::texreg(
    l = list(M1_rob, M2_rob, M3_rob), float.pos = "H",
    digits = 3,
    custom.coef.names = c(
        "(Intercept)" = "Intercept",
        "econ_indic" = "Economic performance",
        "EPI" = "Environmental performance",
        "nr_disasters" = "Number of natural disasters",
        "share_demo" = "Share of demonstrative action"
    ),
    custom.model.names = c("Events", "Participants", "Share of events"),
    omit.coef = "country_iso",
    caption = "Grievances",
    label = "Grievances",
    custom.note = "Robust standard errors in parentheses. %stars",
    custom.gof.rows = list(
        "Country fixed effects" = rep("Yes", 3),
        "Adjusted R sq" = rsqs$adjr2[rsqs$type == "Grievances"],
        "Number of obs." = rep(nrow(reg_data), 3)
    ),
    file = paste0(path, "draft/July_2024/Appendix/Tables/grievances.tex")
)
```


## Opportunities

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
form_events <- as.formula(total_env_protest ~
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_demo +
    country_iso)

form_participants <- as.formula(total_env_part ~
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_demo +
    country_iso)

form_share_events <- as.formula(share_env_protest ~
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_demo +
    country_iso)

M1 <- prais_winsten(
    formula = form_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M1_rob <- coeftest(M1, vcov = vcovPC(M1))

M2 <- prais_winsten(
    formula = form_participants,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M2_rob <- coeftest(M2, vcov = vcovPC(M2))

M3 <- prais_winsten(
    formula = form_share_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M3_rob <- coeftest(M3, vcov = vcovPC(M3))

rsqs$adjr2[rsqs$type == "Opportunities" & rsqs$model == "Events"] <- summary(M1)$adj.r.squared
rsqs$adjr2[rsqs$type == "Opportunities" & rsqs$model == "Participants"] <- summary(M2)$adj.r.squared
rsqs$adjr2[rsqs$type == "Opportunities" & rsqs$model == "Share events"] <- summary(M3)$adj.r.squared

texreg::texreg(
    l = list(M1_rob, M2_rob, M3_rob), float.pos = "H",
    digits = 3,
    custom.coef.names = c(
        "(Intercept)" = "Intercept",
        "gov_env_cmp" = "Gov. position on env. protection",
        "env_prot_sys" = "Parl. position on env. protection",
        "green_in_gov" = "Green party in gov.",
        "green_in_parl" = "Green party in parl.",
        "share_demo" = "Share of demonstrative action"
    ),
    custom.model.names = c("Events", "Participants", "Share of events"),
    omit.coef = "country_iso",
    caption = "Opportunities",
    label = "Opportunities",
    custom.note = "Robust standard errors in parentheses. %stars",
    custom.gof.rows = list(
        "Country fixed effects" = rep("Yes", 3),
        "Adjusted R sq" = rsqs$adjr2[rsqs$type == "Opportunities"],
        "Number of obs." = rep(nrow(reg_data), 3)
    ),
    file = paste0(path, "draft/July_2024/Appendix/Tables/opportunities.tex")
)
```


## Resources

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
form_events <- as.formula(total_env_protest ~
    share_organized +
    NGO_IUCN +
    share_demo +
    country_iso)

form_participants <- as.formula(total_env_part ~
    share_organized +
    NGO_IUCN +
    share_demo +
    country_iso)

form_share_events <- as.formula(share_env_protest ~
    share_organized +
    NGO_IUCN +
    share_demo +
    country_iso)

M1 <- prais_winsten(
    formula = form_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M1_rob <- coeftest(M1, vcov = vcovPC(M1))

M2 <- prais_winsten(
    formula = form_participants,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M2_rob <- coeftest(M2, vcov = vcovPC(M2))

M3 <- prais_winsten(
    formula = form_share_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M3_rob <- coeftest(M3, vcov = vcovPC(M3))

rsqs$adjr2[rsqs$type == "Resources" & rsqs$model == "Events"] <- summary(M1)$adj.r.squared
rsqs$adjr2[rsqs$type == "Resources" & rsqs$model == "Participants"] <- summary(M2)$adj.r.squared
rsqs$adjr2[rsqs$type == "Resources" & rsqs$model == "Share events"] <- summary(M3)$adj.r.squared

texreg(
    l = list(M1_rob, M2_rob, M3_rob), float.pos = "H",
    digits = 3,
    custom.coef.names = c(
        "(Intercept)" = "Intercept",
        "share_organized" = "Share of organized protests",
        "NGO_IUCN" = "Number of env. NGOs",
        "share_demo" = "Share of demonstrative action"
    ),
    custom.model.names = c("Events", "Participants", "Share of events"),
    omit.coef = "country_iso",
    caption = "Resources",
    label = "Resources",
    custom.note = "Robust standard errors in parentheses. %stars",
    custom.gof.rows = list(
        "Country fixed effects" = rep("Yes", 3),
        "Adjusted R sq" = rsqs$adjr2[rsqs$type == "Opportunities"],
        "Number of obs." = rep(nrow(reg_data), 3)
    ),
    file = paste0(path, "draft/July_2024/Appendix/Tables/resources.tex")
)
```


# Relative Importance of grievances, opportunities and resources

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
rsqs$model <- recode(rsqs$model, "Share events" = "Share of events")
rsqs$type <- factor(rsqs$type, levels = c("Grievances", "Opportunities", "Resources"))
rsqs$adjr2 <- as.numeric(rsqs$adjr2)


ggplot(rsqs, aes(x = type, y = adjr2)) +
    geom_bar(stat = "identity", color = "black") +
    facet_grid(~model) +
    labs(x = "Model", y = "Adj. R-square") +
    theme_bw()

ggsave(
    path = paste0(path, "draft/July_2024/Figures/"),
    filename = "Pseudo_R.eps",
    scale = 1.3,
    width = 6, height = 4,
    units = "in", dpi = 300
)
```

# EPI disaggregated

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
form_events <- as.formula(total_env_protest ~
    econ_indic +
    climate_change_index +
    env_health_index + 
    eco_vitality_index +
    nr_disasters +
    share_demo +
    country_iso)

form_participants <- as.formula(total_env_part ~
    econ_indic +
    climate_change_index +
    env_health_index + 
    eco_vitality_index +
    nr_disasters +
    share_demo +
    country_iso)

form_share_events <- as.formula(share_env_protest ~
    econ_indic +
    climate_change_index +
    env_health_index + 
    eco_vitality_index +
    nr_disasters +
    share_demo +
    country_iso)


M1 <- prais_winsten(
    formula = form_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)

# M1_se <- sqrt(diag(vcovPC(M1)))
# coefci(M1, level = 0.95, vcov. = vcovPC(M1))

M1_rob <- coeftest(M1, vcov = vcovPC(M1))

M2 <- prais_winsten(
    formula = form_participants,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)

M2_rob <- coeftest(M2, vcov = vcovPC(M2))

M3 <- prais_winsten(
    formula = form_share_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)

M3_rob <- coeftest(M3, vcov = vcovPC(M3))

rsqs$adjr2[rsqs$type == "Grievances" & rsqs$model == "Events"] <- summary(M1)$adj.r.squared
rsqs$adjr2[rsqs$type == "Grievances" & rsqs$model == "Participants"] <- summary(M2)$adj.r.squared
rsqs$adjr2[rsqs$type == "Grievances" & rsqs$model == "Share events"] <- summary(M3)$adj.r.squared

# Print regression table
texreg::texreg(
    l = list(M1_rob, M2_rob, M3_rob), float.pos = "H",
    digits = 3,
    custom.coef.names = c(
        "(Intercept)" = "Intercept",
        "econ_indic" = "Economic performance",
        "climate_change_index" = "Climate change",
        "env_health_index" = "Environmental health",
        "eco_vitality_index" = "Ecosystem vitality",
        "nr_disasters" = "Number of natural disasters",
        "share_demo" = "Share of demonstrative action"
    ),
    custom.model.names = c("Events", "Participants", "Share of events"),
    omit.coef = "country_iso",
    caption = "Grievances",
    label = "Grievances",
    custom.note = "Robust standard errors in parentheses. %stars",
    custom.gof.rows = list(
        "Country fixed effects" = rep("Yes", 3),
        "Adjusted R sq" = rsqs$adjr2[rsqs$type == "Grievances"],
        "Number of obs." = rep(nrow(reg_data), 3)
    ),
    file = paste0(path, "draft/July_2024/Appendix/Tables/grievances2.tex")
)
```

# Relative Importance of grievances, opportunities and resources

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
rsqs$model <- recode(rsqs$model, "Share events" = "Share of events")
rsqs$type <- factor(rsqs$type, levels = c("Grievances", "Opportunities", "Resources"))
rsqs$adjr2 <- as.numeric(rsqs$adjr2)


ggplot(rsqs, aes(x = type, y = adjr2)) +
    geom_bar(stat = "identity", color = "black") +
    facet_grid(~model) +
    labs(x = "Model", y = "Adj. R-square") +
    theme_bw()

ggsave(
    path = paste0(path, "draft/July_2024/Appendix/Figures/"),
    filename = "Pseudo_R_EPI_disagg.png",
    scale = 1.3,
    width = 6, height = 4,
    units = "in", dpi = 300
)
```

# Regression model: altogether grievances, opportunities and resources

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
form_events <- as.formula(total_env_protest ~
    econ_indic +
    EPI +
    nr_disasters +
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_organized +
    NGO_IUCN +
    share_demo +
    country_iso)

form_participants <- as.formula(total_env_part ~
    econ_indic +
    EPI +
    nr_disasters +
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_organized +
    NGO_IUCN +
    share_demo +
    country_iso)

form_share_events <- as.formula(share_env_protest ~
    econ_indic +
    EPI +
    nr_disasters +
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_organized +
    NGO_IUCN +
    share_demo +
    country_iso)

M1 <- prais_winsten(
    formula = form_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M1_rob <- coeftest(M1, vcov = vcovPC(M1))

M2 <- prais_winsten(
    formula = form_participants,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M2_rob <- coeftest(M2, vcov = vcovPC(M2))

M3 <- prais_winsten(
    formula = form_share_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M3_rob <- coeftest(M3, vcov = vcovPC(M3))

# Print regression results
texreg(
    l = list(M1_rob, M2_rob, M3_rob), float.pos = "H",
    digits = 3,
    custom.coef.names = c(
        "(Intercept)" = "Intercept",
        "econ_indic" = "Economic performance",
        "EPI" = "Environmental performance",
        "nr_disasters" = "Number of natural disasters",
        "gov_env_cmp" = "Gov. position on env. protection",
        "env_prot_sys" = "Parl. position on env. protection",
        "green_in_gov" = "Green party in gov.",
        "green_in_parl" = "Green party in parl.",
        "share_organized" = "Share of organized protests",
        "NGO_IUCN" = "Number of env. NGOs",
        "share_demo" = "Share of demonstrative action"
    ),
    custom.model.names = c("Events", "Participants", "Share of events"),
    omit.coef = "country_iso",
    groups = list(
        "Grievances" = 2:4,
        "Opportunities" = 5:8,
        "Resources" = 9:10
    ),
    caption = "Prais-Winsten regression models of environmental mobilization",
    include.ci = FALSE,
    label = "",
    custom.note = "Robust standard errors in parentheses. %stars",
    custom.gof.rows = list(
        "Country fixed effects" = rep("Yes", 3),
        "Adjusted R sq" = c(summary(M1)$adj.r.squared, summary(M2)$adj.r.squared, summary(M3)$adj.r.squared),
        "Number of obs." = rep(nrow(reg_data), 3)
    ),
    file = paste0(path, "draft/July_2024/Tables/reg_model_full.tex")
)
# "Lagged dependent variables" = c("YES", "YES", "YES"),
```


```{r, echo = F, message = F, warning = F, fig.height = 8, fig.width = 6, fig.align = "center"}
clean_reg_result <- function(model, DV_name) {
    mod_ci <- cbind(
        coefci(model, level = 0.95, vcov. = vcovPC(model)),
        coefci(model, level = 0.84, vcov. = vcovPC(model)),
        coef(model)
    )

    mod_ci <- as.data.frame(mod_ci)
    mod_ci$term <- rownames(mod_ci)
    colnames(mod_ci)[5] <- "estimate"
    mod_ci$dv_name <- DV_name
    rownames(mod_ci) <- NULL

    return(mod_ci)
}

M1_ci <- clean_reg_result(M1, "Events")
M2_ci <- clean_reg_result(M2, "Participants")
M3_ci <- clean_reg_result(M3, "Share of events")
reg_res <- rbind(M1_ci, M2_ci, M3_ci)

reg_res <- reg_res %>%
    filter(!(grepl("country_iso", term))) %>%
    filter(!(grepl("Intercept", term))) %>%
    filter(!term == "share_demo") %>%
    mutate(term = case_when(
        term == "econ_indic" ~ "Econ. perf.",
        term == "EPI" ~ "Environ. perf.",
        term == "nr_disasters" ~ "Natural disasters",
        term == "gov_env_cmp" ~ "Gov. on env. prot.",
        term == "env_prot_sys" ~ "Parl. on env. prot.",
        term == "green_in_gov" ~ "Green prt. in gov.",
        term == "green_in_parl" ~ "Green prt. in parl.",
        term == "share_organized" ~ "Org. protests",
        term == "NGO_IUCN" ~ "Env. NGOs",
        TRUE ~ term,
    )) %>%
    mutate(term = factor(term, levels = c(
        "Econ. perf.", "Environ. perf.",
        "Natural disasters",
        "Gov. on env. prot.", "Parl. on env. prot.",
        "Green prt. in gov.", "Green prt. in parl.",
        "Org. protests", "Env. NGOs"
    ))) %>%
    mutate(term = fct_rev(term)) %>%
    mutate(type = case_when(
        term %in% c(
            "Econ. perf.", 
            "Environ. perf.", "Natural disasters"
        ) ~ "Grievances",
        term %in% c(
            "Gov. on env. prot.", "Parl. on env. prot.",
            "Green prt. in gov.", "Green prt. in parl."
        ) ~ "Opportunities",
        term %in% c("Org. protests", "Env. NGOs") ~ "Resources"
    )) %>%
    mutate(dv_name = fct_rev(dv_name))

space_between_bars <- 0.65

ggplot(reg_res, aes(y = term)) +
    geom_vline(xintercept = 0, linetype = "dashed", color = "gray60") +
    geom_point(aes(x = estimate, shape = dv_name, color = dv_name), position = position_dodge(width = space_between_bars), size = 2) +
    geom_linerange(aes(
        xmin = `2.5 %`, xmax = `97.5 %`, color = dv_name,
        group = dv_name
    ), size = 0.5, position = position_dodge(width = space_between_bars), key_glyph = "path") +
    geom_linerange(aes(
        xmin = `8 %`, xmax = `92 %`, color = dv_name,
        group = dv_name
    ), size = 1.1, position = position_dodge(width = space_between_bars), key_glyph = "path") +
    facet_wrap(~type, scale = "free_y", ncol = 2) +
    xlab("Prais-Winsten Estimate") +
    scale_color_brewer(name = "", palette = "Set2", guide = guide_legend(reverse = TRUE, nrow = 1, byrow = TRUE)) +
    scale_shape_discrete(name = "", guide = guide_legend(reverse = TRUE, nrow = 2, byrow = TRUE)) +
    theme_linedraw() +
    theme(
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.key.width = unit(1, "cm"),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank()
    )

ggsave(
    plot = last_plot(),
    filename = "coefplot_overtime.eps",
    path = paste0(path, "draft/July_2024/Figures/"),
    width = 7,
    height = 6,
    scale = 0.9,
    dpi = 400
)
```

# Regression model same, but EPI disaggregates

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
form_events <- as.formula(total_env_protest ~
    econ_indic +
    climate_change_index +
    env_health_index + 
    eco_vitality_index +
    nr_disasters +
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_organized +
    NGO_IUCN +
    share_demo +
    country_iso)

form_participants <- as.formula(total_env_part ~
    econ_indic +
    climate_change_index +
    env_health_index + 
    eco_vitality_index +
    nr_disasters +
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_organized +
    NGO_IUCN +
    share_demo +
    country_iso)

form_share_events <- as.formula(share_env_protest ~
    econ_indic +
    climate_change_index +
    env_health_index + 
    eco_vitality_index +
    nr_disasters +
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_organized +
    NGO_IUCN +
    share_demo +
    country_iso)

M1 <- prais_winsten(
    formula = form_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M1_rob <- coeftest(M1, vcov = vcovPC(M1))

M2 <- prais_winsten(
    formula = form_participants,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M2_rob <- coeftest(M2, vcov = vcovPC(M2))

M3 <- prais_winsten(
    formula = form_share_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M3_rob <- coeftest(M3, vcov = vcovPC(M3))

# Print regression results
texreg(
    l = list(M1_rob, M2_rob, M3_rob), float.pos = "H",
    digits = 3,
    custom.coef.names = c(
        "(Intercept)" = "Intercept",
        "econ_indic" = "Economic performance",
        "climate_change_index" = "Climate change",
        "env_health_index" = "Environmental health",
        "eco_vitality_index" = "Ecosystem vitality",
        "nr_disasters" = "Number of natural disasters",
        "gov_env_cmp" = "Gov. position on env. protection",
        "env_prot_sys" = "Parl. position on env. protection",
        "green_in_gov" = "Green party in gov.",
        "green_in_parl" = "Green party in parl.",
        "share_organized" = "Share of organized protests",
        "NGO_IUCN" = "Number of env. NGOs",
        "share_demo" = "Share of demonstrative action"
    ),
    custom.model.names = c("Events", "Participants", "Share of events"),
    omit.coef = "country_iso",
    groups = list(
        "Grievances" = 2:6,
        "Opportunities" = 7:10,
        "Resources" = 11:12
    ),
    caption = "Prais-Winsten regression models of environmental mobilization",
    include.ci = FALSE,
    label = "",
    custom.note = "Robust standard errors in parentheses. %stars",
    custom.gof.rows = list(
        "Country fixed effects" = rep("Yes", 3),
        "Adjusted R sq" = c(summary(M1)$adj.r.squared, summary(M2)$adj.r.squared, summary(M3)$adj.r.squared),
        "Number of obs." = rep(nrow(reg_data), 3)
    ),
    file = paste0(path, "draft/July_2024/Appendix/Tables/reg_model_full2.tex")
)
# "Lagged dependent variables" = c("YES", "YES", "YES"),
```


```{r, echo = F, message = F, warning = F, fig.height = 8, fig.width = 6, fig.align = "center"}
clean_reg_result <- function(model, DV_name) {
    mod_ci <- cbind(
        coefci(model, level = 0.95, vcov. = vcovPC(model)),
        coefci(model, level = 0.84, vcov. = vcovPC(model)),
        coef(model)
    )

    mod_ci <- as.data.frame(mod_ci)
    mod_ci$term <- rownames(mod_ci)
    colnames(mod_ci)[5] <- "estimate"
    mod_ci$dv_name <- DV_name
    rownames(mod_ci) <- NULL

    return(mod_ci)
}

M1_ci <- clean_reg_result(M1, "Events")
M2_ci <- clean_reg_result(M2, "Participants")
M3_ci <- clean_reg_result(M3, "Share of events")
reg_res <- rbind(M1_ci, M2_ci, M3_ci)

reg_res <- reg_res %>%
    filter(!(grepl("country_iso", term))) %>%
    filter(!(grepl("Intercept", term))) %>%
    filter(!term == "share_demo") %>%
    mutate(term = case_when(
        term == "econ_indic" ~ "Econ. perf.",
        term == "EPI" ~ "Environ. perf.",
        term == "climate_change_index" ~ "Climate change",
        term == "env_health_index" ~ "Environmental health",
        term == "eco_vitality_index" ~ "Ecosystem vitality",
        term == "nr_disasters" ~ "Natural disasters",
        term == "gov_env_cmp" ~ "Gov. on env. prot.",
        term == "env_prot_sys" ~ "Parl. on env. prot.",
        term == "green_in_gov" ~ "Green prt. in gov.",
        term == "green_in_parl" ~ "Green prt. in parl.",
        term == "share_organized" ~ "Org. protests",
        term == "NGO_IUCN" ~ "Env. NGOs",
        TRUE ~ term,
    )) %>%
    mutate(term = factor(term, levels = c(
        "Econ. perf.", "Climate change", 
        "Environmental health", "Ecosystem vitality",
        "Natural disasters",
        "Gov. on env. prot.", "Parl. on env. prot.",
        "Green prt. in gov.", "Green prt. in parl.",
        "Org. protests", "Env. NGOs"
    ))) %>%
    mutate(term = fct_rev(term)) %>%
    mutate(type = case_when(
        term %in% c(
            "Econ. perf.", 
            "Climate change", "Environmental health",
            "Ecosystem vitality", "Natural disasters"
        ) ~ "Grievances",
        term %in% c(
            "Gov. on env. prot.", "Parl. on env. prot.",
            "Green prt. in gov.", "Green prt. in parl."
        ) ~ "Opportunities",
        term %in% c("Org. protests", "Env. NGOs") ~ "Resources"
    )) %>%
    mutate(dv_name = fct_rev(dv_name))

space_between_bars <- 0.65

ggplot(reg_res, aes(y = term)) +
    geom_vline(xintercept = 0, linetype = "dashed", color = "gray60") +
    geom_point(aes(x = estimate, shape = dv_name, color = dv_name), position = position_dodge(width = space_between_bars), size = 2) +
    geom_linerange(aes(
        xmin = `2.5 %`, xmax = `97.5 %`, color = dv_name,
        group = dv_name
    ), size = 0.5, position = position_dodge(width = space_between_bars), key_glyph = "path") +
    geom_linerange(aes(
        xmin = `8 %`, xmax = `92 %`, color = dv_name,
        group = dv_name
    ), size = 1.1, position = position_dodge(width = space_between_bars), key_glyph = "path") +
    facet_wrap(~type, scale = "free_y", ncol = 2) +
    xlab("Prais-Winsten Estimate") +
    scale_color_brewer(name = "", palette = "Set2", guide = guide_legend(reverse = TRUE, nrow = 1, byrow = TRUE)) +
    scale_shape_discrete(name = "", guide = guide_legend(reverse = TRUE, nrow = 2, byrow = TRUE)) +
    theme_linedraw() +
    theme(
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.key.width = unit(1, "cm"),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank()
    )

ggsave(
    plot = last_plot(),
    filename = "coefplot_overtime2.png",
    path = paste0(path, "draft/July_2024/Appendix/Figures/"),
    width = 7,
    height = 6,
    scale = 0.9,
    dpi = 400
)
```



# Regression models - country means

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
reg_data <- subset_reg_data(threshold = 20)

table(is.na(reg_data$year))

dput(sort(countrycode::countrycode(unique(PEA_yearly$country_iso), origin = "iso2c", destination = "country.name")))

reg_data <- reg_data %>%
    select(
        total_env_protest, total_env_part, share_env_protest,
        econ_indic, EPI, climate_change_index, env_health_index, 
        eco_vitality_index, nr_disasters, share_demo,
        gov_env_cmp, env_prot_sys, green_in_gov, green_in_parl,
        share_organized, NGO_IUCN,
        country_iso
    ) %>%
    group_by(country_iso) %>%
    mutate_at(vars(
        total_env_protest, total_env_part, share_env_protest, econ_indic, 
        EPI, climate_change_index, env_health_index, 
        eco_vitality_index, nr_disasters, share_demo,
        gov_env_cmp, env_prot_sys, green_in_gov, green_in_parl,
        share_organized, NGO_IUCN
    ), ~ mean(., na.rm = TRUE)) %>%
    ungroup() %>%
    unique() %>%
    na.omit() %>%
    mutate_at(vars(
        total_env_protest, total_env_part, share_env_protest, econ_indic,
        EPI, climate_change_index, env_health_index, 
        eco_vitality_index, nr_disasters, share_demo,
        gov_env_cmp, env_prot_sys, green_in_gov, green_in_parl,
        share_organized, NGO_IUCN
    ), ~ scale(.))

```

## bivariate model - EPI

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
bivariate_reg <- function(DV, IV) {
    form <- as.formula(paste0(DV, " ~ ", IV))

    model <- lm(form, data = reg_data)
    tidy(model, conf.level = 0.84, conf.int = TRUE)

    mod_ci <- cbind(
        coefci(model, level = 0.95),
        coefci(model, level = 0.84),
        coef(model)
    )

    mod_ci <- as.data.frame(mod_ci)
    mod_ci <- mod_ci[2, ]
    mod_ci$term <- rownames(mod_ci)
    colnames(mod_ci)[5] <- "estimate"
    mod_ci$dv_name <- DV
    rownames(mod_ci) <- NULL

    return(mod_ci)
}

reg_res <- bind_rows(
    bivariate_reg("total_env_protest", "econ_indic"),
    bivariate_reg("total_env_protest", "EPI"),
    bivariate_reg("total_env_protest", "nr_disasters"),
    bivariate_reg("total_env_protest", "gov_env_cmp"),
    bivariate_reg("total_env_protest", "env_prot_sys"),
    bivariate_reg("total_env_protest", "green_in_gov"),
    bivariate_reg("total_env_protest", "green_in_parl"),
    bivariate_reg("total_env_protest", "share_organized"),
    bivariate_reg("total_env_protest", "NGO_IUCN"),
    bivariate_reg("total_env_part", "econ_indic"),
    bivariate_reg("total_env_part", "EPI"),
    bivariate_reg("total_env_part", "nr_disasters"),
    bivariate_reg("total_env_part", "gov_env_cmp"),
    bivariate_reg("total_env_part", "env_prot_sys"),
    bivariate_reg("total_env_part", "green_in_gov"),
    bivariate_reg("total_env_part", "green_in_parl"),
    bivariate_reg("total_env_part", "share_organized"),
    bivariate_reg("total_env_part", "NGO_IUCN"),
    bivariate_reg("share_env_protest", "econ_indic"),
    bivariate_reg("share_env_protest", "EPI"),
    bivariate_reg("share_env_protest", "nr_disasters"),
    bivariate_reg("share_env_protest", "gov_env_cmp"),
    bivariate_reg("share_env_protest", "env_prot_sys"),
    bivariate_reg("share_env_protest", "green_in_gov"),
    bivariate_reg("share_env_protest", "green_in_parl"),
    bivariate_reg("share_env_protest", "share_organized"),
    bivariate_reg("share_env_protest", "NGO_IUCN")
)

reg_res <- reg_res %>%
    mutate(term = case_when(
        term == "econ_indic" ~ "Econ. perf.",
        term == "EPI" ~ "Environ. perf.",
        term == "nr_disasters" ~ "Natural disasters",
        term == "gov_env_cmp" ~ "Gov. on env. prot.",
        term == "env_prot_sys" ~ "Parl. on env. prot.",
        term == "green_in_gov" ~ "Green prt. in gov.",
        term == "green_in_parl" ~ "Green prt. in parl.",
        term == "share_organized" ~ "Org. protests",
        term == "NGO_IUCN" ~ "Env. NGOs",
        TRUE ~ term,
    )) %>%
    mutate(term = factor(term, levels = c(
        "Econ. perf.", 
        "Environ. perf.", "Natural disasters",
        "Gov. on env. prot.", "Parl. on env. prot.",
        "Green prt. in gov.", "Green prt. in parl.",
        "Org. protests", "Env. NGOs"
    ))) %>%
    mutate(term = fct_rev(term)) %>%
    mutate(type = case_when(
        term %in% c(
            "Econ. perf.", 
            "Environ. perf.", "Natural disasters"
        ) ~ "Grievances",
        term %in% c(
            "Gov. on env. prot.", "Parl. on env. prot.",
            "Green prt. in gov.", "Green prt. in parl."
        ) ~ "Opportunities",
        term %in% c("Org. protests", "Env. NGOs") ~ "Resources"
    )) %>%
    mutate(dv_name = case_when(
        dv_name == "total_env_protest" ~ "Events",
        dv_name == "total_env_part" ~ "Participants",
        dv_name == "share_env_protest" ~ "Share of events"
    )) %>%
    mutate(dv_name = fct_rev(dv_name))

space_between_bars <- 0.65

ggplot(reg_res, aes(y = term)) +
    geom_vline(xintercept = 0, linetype = "dashed", color = "gray60") +
    geom_point(aes(x = estimate, shape = dv_name, color = dv_name), position = position_dodge(width = space_between_bars), size = 2) +
    geom_linerange(aes(
        xmin = `2.5 %`, xmax = `97.5 %`, color = dv_name,
        group = dv_name
    ), size = 0.5, position = position_dodge(width = space_between_bars), key_glyph = "path") +
    geom_linerange(aes(
        xmin = `8 %`, xmax = `92 %`, color = dv_name,
        group = dv_name
    ), size = 1.1, position = position_dodge(width = space_between_bars), key_glyph = "path") +
    facet_wrap(~type, scale = "free_y", ncol = 2) +
    xlab("OLS Estimates") +
    scale_color_brewer(name = "Dependent variables:", palette = "Set2", guide = guide_legend(reverse = TRUE, nrow = 1, byrow = TRUE)) +
    scale_shape_discrete(name = "Dependent variables:", guide = guide_legend(reverse = TRUE, nrow = 2, byrow = TRUE)) +
    theme_linedraw() +
    theme(
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.key.width = unit(1, "cm"),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank()
    )

ggsave(
    plot = last_plot(),
    filename = "coefplot_countries_bivariate.png",
    path = paste0(path, "draft/July_2024/Appendix/Figures/"),
    width = 7,
    height = 6,
    scale = 0.9,
    dpi = 400
)
```

## bivariate model - EPI disaggregated

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
bivariate_reg <- function(DV, IV) {
    form <- as.formula(paste0(DV, " ~ ", IV))

    model <- lm(form, data = reg_data)
    tidy(model, conf.level = 0.84, conf.int = TRUE)

    mod_ci <- cbind(
        coefci(model, level = 0.95),
        coefci(model, level = 0.84),
        coef(model)
    )

    mod_ci <- as.data.frame(mod_ci)
    mod_ci <- mod_ci[2, ]
    mod_ci$term <- rownames(mod_ci)
    colnames(mod_ci)[5] <- "estimate"
    mod_ci$dv_name <- DV
    rownames(mod_ci) <- NULL

    return(mod_ci)
}

reg_res <- bind_rows(
    bivariate_reg("total_env_protest", "econ_indic"),
    bivariate_reg("total_env_protest", "climate_change_index"),
    bivariate_reg("total_env_protest", "env_health_index"),
    bivariate_reg("total_env_protest", "eco_vitality_index"),
    bivariate_reg("total_env_protest", "nr_disasters"),
    bivariate_reg("total_env_protest", "gov_env_cmp"),
    bivariate_reg("total_env_protest", "env_prot_sys"),
    bivariate_reg("total_env_protest", "green_in_gov"),
    bivariate_reg("total_env_protest", "green_in_parl"),
    bivariate_reg("total_env_protest", "share_organized"),
    bivariate_reg("total_env_protest", "NGO_IUCN"),
    bivariate_reg("total_env_part", "econ_indic"),
    bivariate_reg("total_env_part", "climate_change_index"),
    bivariate_reg("total_env_part", "env_health_index"),
    bivariate_reg("total_env_part", "eco_vitality_index"),
    bivariate_reg("total_env_part", "nr_disasters"),
    bivariate_reg("total_env_part", "gov_env_cmp"),
    bivariate_reg("total_env_part", "env_prot_sys"),
    bivariate_reg("total_env_part", "green_in_gov"),
    bivariate_reg("total_env_part", "green_in_parl"),
    bivariate_reg("total_env_part", "share_organized"),
    bivariate_reg("total_env_part", "NGO_IUCN"),
    bivariate_reg("share_env_protest", "econ_indic"),
    bivariate_reg("share_env_protest", "climate_change_index"),
    bivariate_reg("share_env_protest", "env_health_index"),
    bivariate_reg("share_env_protest", "eco_vitality_index"),
    bivariate_reg("share_env_protest", "nr_disasters"),
    bivariate_reg("share_env_protest", "gov_env_cmp"),
    bivariate_reg("share_env_protest", "env_prot_sys"),
    bivariate_reg("share_env_protest", "green_in_gov"),
    bivariate_reg("share_env_protest", "green_in_parl"),
    bivariate_reg("share_env_protest", "share_organized"),
    bivariate_reg("share_env_protest", "NGO_IUCN")
)

reg_res <- reg_res %>%
    mutate(term = case_when(
        term == "econ_indic" ~ "Econ. perf.",
        term == "climate_change_index" ~ "Climate change",
        term == "env_health_index" ~ "Environmental health",
        term == "eco_vitality_index" ~ "Ecosystem vitality",
        term == "nr_disasters" ~ "Natural disasters",
        term == "gov_env_cmp" ~ "Gov. on env. prot.",
        term == "env_prot_sys" ~ "Parl. on env. prot.",
        term == "green_in_gov" ~ "Green prt. in gov.",
        term == "green_in_parl" ~ "Green prt. in parl.",
        term == "share_organized" ~ "Org. protests",
        term == "NGO_IUCN" ~ "Env. NGOs",
        TRUE ~ term,
    )) %>%
    mutate(term = factor(term, levels = c(
        "Econ. perf.", 
        "Climate change", "Environmental health",
        "Ecosystem vitality", "Natural disasters",
        "Gov. on env. prot.", "Parl. on env. prot.",
        "Green prt. in gov.", "Green prt. in parl.",
        "Org. protests", "Env. NGOs"
    ))) %>%
    mutate(term = fct_rev(term)) %>%
    mutate(type = case_when(
        term %in% c(
            "Econ. perf.", 
            "Climate change", "Environmental health",
            "Ecosystem vitality", "Natural disasters"
        ) ~ "Grievances",
        term %in% c(
            "Gov. on env. prot.", "Parl. on env. prot.",
            "Green prt. in gov.", "Green prt. in parl."
        ) ~ "Opportunities",
        term %in% c("Org. protests", "Env. NGOs") ~ "Resources"
    )) %>%
    mutate(dv_name = case_when(
        dv_name == "total_env_protest" ~ "Events",
        dv_name == "total_env_part" ~ "Participants",
        dv_name == "share_env_protest" ~ "Share of events"
    )) %>%
    mutate(dv_name = fct_rev(dv_name))

space_between_bars <- 0.65

ggplot(reg_res, aes(y = term)) +
    geom_vline(xintercept = 0, linetype = "dashed", color = "gray60") +
    geom_point(aes(x = estimate, shape = dv_name, color = dv_name), position = position_dodge(width = space_between_bars), size = 2) +
    geom_linerange(aes(
        xmin = `2.5 %`, xmax = `97.5 %`, color = dv_name,
        group = dv_name
    ), size = 0.5, position = position_dodge(width = space_between_bars), key_glyph = "path") +
    geom_linerange(aes(
        xmin = `8 %`, xmax = `92 %`, color = dv_name,
        group = dv_name
    ), size = 1.1, position = position_dodge(width = space_between_bars), key_glyph = "path") +
    facet_wrap(~type, scale = "free_y", ncol = 2) +
    xlab("OLS Estimates") +
    scale_color_brewer(name = "Dependent variables:", palette = "Set2", guide = guide_legend(reverse = TRUE, nrow = 1, byrow = TRUE)) +
    scale_shape_discrete(name = "Dependent variables:", guide = guide_legend(reverse = TRUE, nrow = 2, byrow = TRUE)) +
    theme_linedraw() +
    theme(
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.key.width = unit(1, "cm"),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank()
    )

ggsave(
    plot = last_plot(),
    filename = "coefplot_countries_bivariate2.png",
    path = paste0(path, "draft/July_2024/Appendix/Figures/"),
    width = 7,
    height = 6,
    scale = 0.9,
    dpi = 400
)
```

# multivariate model

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
form_events <- as.formula(total_env_protest ~
    econ_indic +
    EPI +
    nr_disasters +
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_organized +
    NGO_IUCN)
# +
# share_demo

form_participants <- as.formula(total_env_part ~
    econ_indic +
    EPI +
    nr_disasters +
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_organized +
    NGO_IUCN)

form_share_events <- as.formula(share_env_protest ~
    econ_indic +
    EPI +
    nr_disasters +
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_organized +
    NGO_IUCN)

M1 <- lm(
    formula = form_events,
    data = reg_data
)

M2 <- lm(
    formula = form_participants,
    data = reg_data
)

M3 <- lm(
    formula = form_share_events,
    data = reg_data,
)

# Print regression results
texreg(
    l = list(M1, M2, M3), float.pos = "H",
    digits = 3,
    custom.coef.names = c(
        "(Intercept)" = "Intercept",
        "econ_indic" = "Economic performance",
        "EPI" = "Environmental performance",
        "nr_disasters" = "Number of natural disasters",
        "gov_env_cmp" = "Gov. position on env. protection",
        "env_prot_sys" = "Parl. position on env. protection",
        "green_in_gov" = "Green party in gov.",
        "green_in_parl" = "Green party in parl.",
        "share_organized" = "Share of organized protests",
        "NGO_IUCN" = "Number of env. NGOs"
        # "share_demo" = "Share of demonstrative action"
    ),
    custom.model.names = c("Events", "Participants", "Share of events"),
    omit.coef = "country_iso",
    groups = list(
        "Grievances" = 2:4,
        "Opportunities" = 5:8,
        "Resources" = 9:10
    ),
    caption = "Country-level OLS multivariate regression model",
    label = "",
    file = paste0(path, "draft/July_2024/Appendix/Tables/reg_model_cntry_multivar.tex")
)
```


```{r, echo = F, message = F, warning = F, fig.height = 8, fig.width = 6, fig.align = "center"}
clean_reg_result <- function(model, DV_name) {
    mod_ci <- cbind(
        coefci(model, level = 0.95),
        coefci(model, level = 0.84),
        coef(model)
    )

    mod_ci <- as.data.frame(mod_ci)
    mod_ci$term <- rownames(mod_ci)
    colnames(mod_ci)[5] <- "estimate"
    mod_ci$dv_name <- DV_name
    rownames(mod_ci) <- NULL

    return(mod_ci)
}

M1_ci <- clean_reg_result(M1, "Events")
M2_ci <- clean_reg_result(M2, "Participants")
M3_ci <- clean_reg_result(M3, "Share of events")
reg_res <- rbind(M1_ci, M2_ci, M3_ci)

reg_res <- reg_res %>%
    filter(!(grepl("country_iso", term))) %>%
    filter(!(grepl("Intercept", term))) %>%
    filter(!term == "share_demo") %>%
    mutate(term = case_when(
        term == "econ_indic" ~ "Econ. perf.",
        term == "EPI" ~ "Environ. perf.",
        term == "nr_disasters" ~ "Natural disasters",
        term == "gov_env_cmp" ~ "Gov. on env. prot.",
        term == "env_prot_sys" ~ "Parl. on env. prot.",
        term == "green_in_gov" ~ "Green prt. in gov.",
        term == "green_in_parl" ~ "Green prt. in parl.",
        term == "share_organized" ~ "Org. protests",
        term == "NGO_IUCN" ~ "Env. NGOs",
        TRUE ~ term,
    )) %>%
    mutate(term = factor(term, levels = c(
        "Econ. perf.", 
        "Environ. perf.", "Natural disasters",
        "Gov. on env. prot.", "Parl. on env. prot.",
        "Green prt. in gov.", "Green prt. in parl.",
        "Org. protests", "Env. NGOs"
    ))) %>%
    mutate(term = fct_rev(term)) %>%
    mutate(type = case_when(
        term %in% c(
            "Econ. perf.", 
            "Environ. perf.", "Natural disasters"
        ) ~ "Grievances",
        term %in% c(
            "Gov. on env. prot.", "Parl. on env. prot.",
            "Green prt. in gov.", "Green prt. in parl."
        ) ~ "Opportunities",
        term %in% c("Org. protests", "Env. NGOs") ~ "Resources"
    )) %>%
    mutate(dv_name = fct_rev(dv_name))

space_between_bars <- 0.65

ggplot(reg_res, aes(y = term)) +
    geom_vline(xintercept = 0, linetype = "dashed", color = "gray60") +
    geom_point(aes(x = estimate, shape = dv_name, color = dv_name), position = position_dodge(width = space_between_bars), size = 2) +
    geom_linerange(aes(
        xmin = `2.5 %`, xmax = `97.5 %`, color = dv_name,
        group = dv_name
    ), size = 0.5, position = position_dodge(width = space_between_bars), key_glyph = "path") +
    geom_linerange(aes(
        xmin = `8 %`, xmax = `92 %`, color = dv_name,
        group = dv_name
    ), size = 1.1, position = position_dodge(width = space_between_bars), key_glyph = "path") +
    facet_wrap(~type, scale = "free_y", ncol = 2) +
    xlab("OLS Estimates") +
    scale_color_brewer(name = "", palette = "Set2", guide = guide_legend(reverse = TRUE, nrow = 1, byrow = TRUE)) +
    scale_shape_discrete(name = "", guide = guide_legend(reverse = TRUE, nrow = 2, byrow = TRUE)) +
    theme_linedraw() +
    theme(
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.key.width = unit(1, "cm"),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank()
    )

ggsave(
    plot = last_plot(),
    filename = "coefplot_countries_multivariate.png",
    path = paste0(path, "draft/July_2024/Appendix/Figures/"),
    width = 7,
    height = 6,
    scale = 0.9,
    dpi = 400
)
```


# multivariate model - EPI disaggregated

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
form_events <- as.formula(total_env_protest ~
    econ_indic +
    climate_change_index +
    env_health_index + 
    eco_vitality_index +
    nr_disasters +
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_organized +
    NGO_IUCN)
# +
# share_demo

form_participants <- as.formula(total_env_part ~
    econ_indic +
    climate_change_index +
    env_health_index + 
    eco_vitality_index +
    nr_disasters +
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_organized +
    NGO_IUCN)

form_share_events <- as.formula(share_env_protest ~
    econ_indic +
    climate_change_index +
    env_health_index + 
    eco_vitality_index +
    nr_disasters +
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_organized +
    NGO_IUCN)

M1 <- lm(
    formula = form_events,
    data = reg_data
)

M2 <- lm(
    formula = form_participants,
    data = reg_data
)

M3 <- lm(
    formula = form_share_events,
    data = reg_data,
)

# Print regression results
texreg(
    l = list(M1, M2, M3), float.pos = "H",
    digits = 3,
    custom.coef.names = c(
        "(Intercept)" = "Intercept",
        "econ_indic" = "Economic performance",
        "climate_change_index" = "Climate change",
        "env_health_index" = "Environmental health",
        "eco_vitality_index" = "Ecosystem vitality",
        "nr_disasters" = "Number of natural disasters",
        "gov_env_cmp" = "Gov. position on env. protection",
        "env_prot_sys" = "Parl. position on env. protection",
        "green_in_gov" = "Green party in gov.",
        "green_in_parl" = "Green party in parl.",
        "share_organized" = "Share of organized protests",
        "NGO_IUCN" = "Number of env. NGOs"
        # "share_demo" = "Share of demonstrative action"
    ),
    custom.model.names = c("Events", "Participants", "Share of events"),
    omit.coef = "country_iso",
    groups = list(
        "Grievances" = 2:6,
        "Opportunities" = 7:10,
        "Resources" = 11:12
    ),
    caption = "Country-level OLS multivariate regression model",
    label = "",
    file = paste0(path, "draft/July_2024/Appendix/Tables/reg_model_cntry_multivar2.tex")
)
```


```{r, echo = F, message = F, warning = F, fig.height = 8, fig.width = 6, fig.align = "center"}
clean_reg_result <- function(model, DV_name) {
    mod_ci <- cbind(
        coefci(model, level = 0.95),
        coefci(model, level = 0.84),
        coef(model)
    )

    mod_ci <- as.data.frame(mod_ci)
    mod_ci$term <- rownames(mod_ci)
    colnames(mod_ci)[5] <- "estimate"
    mod_ci$dv_name <- DV_name
    rownames(mod_ci) <- NULL

    return(mod_ci)
}

M1_ci <- clean_reg_result(M1, "Events")
M2_ci <- clean_reg_result(M2, "Participants")
M3_ci <- clean_reg_result(M3, "Share of events")
reg_res <- rbind(M1_ci, M2_ci, M3_ci)

reg_res <- reg_res %>%
    filter(!(grepl("country_iso", term))) %>%
    filter(!(grepl("Intercept", term))) %>%
    filter(!term == "share_demo") %>%
    mutate(term = case_when(
        term == "econ_indic" ~ "Econ. perf.",
        term == "climate_change_index" ~ "Climate change",
        term == "env_health_index" ~ "Environmental health",
        term == "eco_vitality_index" ~ "Ecosystem vitality",
        term == "nr_disasters" ~ "Natural disasters",
        term == "gov_env_cmp" ~ "Gov. on env. prot.",
        term == "env_prot_sys" ~ "Parl. on env. prot.",
        term == "green_in_gov" ~ "Green prt. in gov.",
        term == "green_in_parl" ~ "Green prt. in parl.",
        term == "share_organized" ~ "Org. protests",
        term == "NGO_IUCN" ~ "Env. NGOs",
        TRUE ~ term,
    )) %>%
    mutate(term = factor(term, levels = c(
        "Econ. perf.", "Climate change",
        "Environmental health", "Ecosystem vitality",
        "Natural disasters",
        "Gov. on env. prot.", "Parl. on env. prot.",
        "Green prt. in gov.", "Green prt. in parl.",
        "Org. protests", "Env. NGOs"
    ))) %>%
    mutate(term = fct_rev(term)) %>%
    mutate(type = case_when(
        term %in% c(
            "Econ. perf.", "Climate change",
            "Environmental health", "Ecosystem vitality", 
            "Natural disasters"
        ) ~ "Grievances",
        term %in% c(
            "Gov. on env. prot.", "Parl. on env. prot.",
            "Green prt. in gov.", "Green prt. in parl."
        ) ~ "Opportunities",
        term %in% c("Org. protests", "Env. NGOs") ~ "Resources"
    )) %>%
    mutate(dv_name = fct_rev(dv_name))

space_between_bars <- 0.65

ggplot(reg_res, aes(y = term)) +
    geom_vline(xintercept = 0, linetype = "dashed", color = "gray60") +
    geom_point(aes(x = estimate, shape = dv_name, color = dv_name), position = position_dodge(width = space_between_bars), size = 2) +
    geom_linerange(aes(
        xmin = `2.5 %`, xmax = `97.5 %`, color = dv_name,
        group = dv_name
    ), size = 0.5, position = position_dodge(width = space_between_bars), key_glyph = "path") +
    geom_linerange(aes(
        xmin = `8 %`, xmax = `92 %`, color = dv_name,
        group = dv_name
    ), size = 1.1, position = position_dodge(width = space_between_bars), key_glyph = "path") +
    facet_wrap(~type, scale = "free_y", ncol = 2) +
    xlab("OLS Estimates") +
    scale_color_brewer(name = "", palette = "Set2", guide = guide_legend(reverse = TRUE, nrow = 1, byrow = TRUE)) +
    scale_shape_discrete(name = "", guide = guide_legend(reverse = TRUE, nrow = 2, byrow = TRUE)) +
    theme_linedraw() +
    theme(
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.key.width = unit(1, "cm"),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank()
    )

ggsave(
    plot = last_plot(),
    filename = "coefplot_countries_multivariate2.png",
    path = paste0(path, "draft/July_2024/Appendix/Figures/"),
    width = 7,
    height = 6,
    scale = 0.9,
    dpi = 400
)
```


# Replicating the Prais Winsten regression on a full sample

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}

reg_data <- subset_reg_data(threshold = 0)

table(is.na(reg_data$year))

dput(sort(countrycode::countrycode(unique(PEA_yearly$country_iso), origin = "iso2c", destination = "country.name")))

reg_data <- reg_data %>%
    select(
        total_env_protest, total_env_part, share_env_protest,
        econ_indic, EPI, nr_disasters, share_demo,
        gov_env_cmp, env_prot_sys, green_in_gov, green_in_parl,
        share_organized, NGO_IUCN,
        country_iso, year
    ) 

Amelia::missmap(reg_data)

reg_data <- reg_data %>%
    na.omit() %>%
    mutate_at(vars(
        total_env_protest, total_env_part, share_env_protest, econ_indic, 
        EPI, nr_disasters, share_demo,
        gov_env_cmp, env_prot_sys, green_in_gov, green_in_parl,
        share_organized, NGO_IUCN
    ), ~ scale(.))

rsqs <- as.data.frame(expand.grid(
    type = c("Grievances", "Opportunities", "Resources"),
    model = c("Events", "Participants", "Share events"),
    adjr2 = as.numeric(NA)
))
```



## Grievances

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
form_events <- as.formula(total_env_protest ~
    econ_indic +
    EPI +
    nr_disasters +
    share_demo +
    country_iso)

form_participants <- as.formula(total_env_part ~
    econ_indic +
    EPI +
    nr_disasters +
    share_demo +
    country_iso)

form_share_events <- as.formula(share_env_protest ~
    econ_indic +
    EPI +
    nr_disasters +
    share_demo +
    country_iso)


M1 <- prais_winsten(
    formula = form_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)

# M1_se <- sqrt(diag(vcovPC(M1)))
# coefci(M1, level = 0.95, vcov. = vcovPC(M1))

M1_rob <- coeftest(M1, vcov = vcovPC(M1))

M2 <- prais_winsten(
    formula = form_participants,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)

M2_rob <- coeftest(M2, vcov = vcovPC(M2))

M3 <- prais_winsten(
    formula = form_share_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)

M3_rob <- coeftest(M3, vcov = vcovPC(M3))

rsqs$adjr2[rsqs$type == "Grievances" & rsqs$model == "Events"] <- summary(M1)$adj.r.squared
rsqs$adjr2[rsqs$type == "Grievances" & rsqs$model == "Participants"] <- summary(M2)$adj.r.squared
rsqs$adjr2[rsqs$type == "Grievances" & rsqs$model == "Share events"] <- summary(M3)$adj.r.squared

# Print regression table
texreg::texreg(
    l = list(M1_rob, M2_rob, M3_rob), float.pos = "H",
    digits = 3,
    custom.coef.names = c(
        "(Intercept)" = "Intercept",
        "econ_indic" = "Economic performance",
        "EPI" = "Environmental performance",
        "nr_disasters" = "Number of natural disasters",
        "share_demo" = "Share of demonstrative action"
    ),
    custom.model.names = c("Events", "Participants", "Share of events"),
    omit.coef = "country_iso",
    caption = "Grievances",
    label = "Grievances",
    custom.note = "Robust standard errors in parentheses. %stars",
    custom.gof.rows = list(
        "Country fixed effects" = rep("Yes", 3),
        "Adjusted R sq" = rsqs$adjr2[rsqs$type == "Grievances"],
        "Number of obs." = rep(nrow(reg_data), 3)
    ),
    file = paste0(path, "draft/July_2024/Appendix/Tables/grievances_full_sample.tex")
)
```

## Opportunities

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
form_events <- as.formula(total_env_protest ~
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_demo +
    country_iso)

form_participants <- as.formula(total_env_part ~
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_demo +
    country_iso)

form_share_events <- as.formula(share_env_protest ~
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_demo +
    country_iso)

M1 <- prais_winsten(
    formula = form_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M1_rob <- coeftest(M1, vcov = vcovPC(M1))

M2 <- prais_winsten(
    formula = form_participants,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M2_rob <- coeftest(M2, vcov = vcovPC(M2))

M3 <- prais_winsten(
    formula = form_share_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M3_rob <- coeftest(M3, vcov = vcovPC(M3))

rsqs$adjr2[rsqs$type == "Opportunities" & rsqs$model == "Events"] <- summary(M1)$adj.r.squared
rsqs$adjr2[rsqs$type == "Opportunities" & rsqs$model == "Participants"] <- summary(M2)$adj.r.squared
rsqs$adjr2[rsqs$type == "Opportunities" & rsqs$model == "Share events"] <- summary(M3)$adj.r.squared

texreg::texreg(
    l = list(M1_rob, M2_rob, M3_rob), float.pos = "H",
    digits = 3,
    custom.coef.names = c(
        "(Intercept)" = "Intercept",
        "gov_env_cmp" = "Gov. position on env. protection",
        "env_prot_sys" = "Parl. position on env. protection",
        "green_in_gov" = "Green party in gov.",
        "green_in_parl" = "Green party in parl.",
        "share_demo" = "Share of demonstrative action"
    ),
    custom.model.names = c("Events", "Participants", "Share of events"),
    omit.coef = "country_iso",
    caption = "Opportunities",
    label = "Opportunities",
    custom.note = "Robust standard errors in parentheses. %stars",
    custom.gof.rows = list(
        "Country fixed effects" = rep("Yes", 3),
        "Adjusted R sq" = rsqs$adjr2[rsqs$type == "Opportunities"],
        "Number of obs." = rep(nrow(reg_data), 3)
    ),
    file = paste0(path, "draft/July_2024/Appendix/Tables/opportunities_full_sample.tex")
)
```


## Resources

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
form_events <- as.formula(total_env_protest ~
    share_organized +
    NGO_IUCN +
    share_demo +
    country_iso)

form_participants <- as.formula(total_env_part ~
    share_organized +
    NGO_IUCN +
    share_demo +
    country_iso)

form_share_events <- as.formula(share_env_protest ~
    share_organized +
    NGO_IUCN +
    share_demo +
    country_iso)

M1 <- prais_winsten(
    formula = form_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M1_rob <- coeftest(M1, vcov = vcovPC(M1))

M2 <- prais_winsten(
    formula = form_participants,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M2_rob <- coeftest(M2, vcov = vcovPC(M2))

M3 <- prais_winsten(
    formula = form_share_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M3_rob <- coeftest(M3, vcov = vcovPC(M3))

rsqs$adjr2[rsqs$type == "Resources" & rsqs$model == "Events"] <- summary(M1)$adj.r.squared
rsqs$adjr2[rsqs$type == "Resources" & rsqs$model == "Participants"] <- summary(M2)$adj.r.squared
rsqs$adjr2[rsqs$type == "Resources" & rsqs$model == "Share events"] <- summary(M3)$adj.r.squared

texreg(
    l = list(M1_rob, M2_rob, M3_rob), float.pos = "H",
    digits = 3,
    custom.coef.names = c(
        "(Intercept)" = "Intercept",
        "share_organized" = "Share of organized protests",
        "NGO_IUCN" = "Number of env. NGOs",
        "share_demo" = "Share of demonstrative action"
    ),
    custom.model.names = c("Events", "Participants", "Share of events"),
    omit.coef = "country_iso",
    caption = "Resources",
    label = "Resources",
    custom.note = "Robust standard errors in parentheses. %stars",
    custom.gof.rows = list(
        "Country fixed effects" = rep("Yes", 3),
        "Adjusted R sq" = rsqs$adjr2[rsqs$type == "Opportunities"],
        "Number of obs." = rep(nrow(reg_data), 3)
    ),
    file = paste0(path, "draft/July_2024/Appendix/Tables/resources_full_sample.tex")
)
```


# Relative Importance of grievances, opportunities and resources

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
rsqs$model <- recode(rsqs$model, "Share events" = "Share of events")
rsqs$type <- factor(rsqs$type, levels = c("Grievances", "Opportunities", "Resources"))
rsqs$adjr2 <- as.numeric(rsqs$adjr2)


ggplot(rsqs, aes(x = type, y = adjr2)) +
    geom_bar(stat = "identity", color = "black") +
    facet_grid(~model) +
    labs(x = "Model", y = "Adj. R-square") +
    theme_bw()

ggsave(
    path = paste0(path, "draft/July_2024/Appendix/Figures/"),
    filename = "Pseudo_R_full_sample.png",
    scale = 1.3,
    width = 6, height = 4,
    units = "in", dpi = 300
)
```

# Regression model: altogether grievances, opportunities and resources

```{r, results = 'asis', include=TRUE, echo=FALSE, warning=FALSE, message=FALSE}
form_events <- as.formula(total_env_protest ~
    econ_indic +
    EPI +
    nr_disasters +
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_organized +
    NGO_IUCN +
    share_demo +
    country_iso)

form_participants <- as.formula(total_env_part ~
    econ_indic +
    EPI +
    nr_disasters +
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_organized +
    NGO_IUCN +
    share_demo +
    country_iso)

form_share_events <- as.formula(share_env_protest ~
    econ_indic +
    EPI +
    nr_disasters +
    gov_env_cmp +
    env_prot_sys +
    green_in_gov +
    green_in_parl +
    share_organized +
    NGO_IUCN +
    share_demo +
    country_iso)

M1 <- prais_winsten(
    formula = form_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M1_rob <- coeftest(M1, vcov = vcovPC(M1))

M2 <- prais_winsten(
    formula = form_participants,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M2_rob <- coeftest(M2, vcov = vcovPC(M2))

M3 <- prais_winsten(
    formula = form_share_events,
    data = reg_data,
    panelwise = TRUE,
    index = c("country_iso", "year")
)
M3_rob <- coeftest(M3, vcov = vcovPC(M3))

# Print regression results
texreg(
    l = list(M1_rob, M2_rob, M3_rob), float.pos = "H",
    digits = 3,
    custom.coef.names = c(
        "(Intercept)" = "Intercept",
        "econ_indic" = "Economic performance",
        "EPI" = "Environmental performance ",
        "nr_disasters" = "Number of natural disasters",
        "gov_env_cmp" = "Gov. position on env. protection",
        "env_prot_sys" = "Parl. position on env. protection",
        "green_in_gov" = "Green party in gov.",
        "green_in_parl" = "Green party in parl.",
        "share_organized" = "Share of organized protests",
        "NGO_IUCN" = "Number of env. NGOs",
        "share_demo" = "Share of demonstrative action"
    ),
    custom.model.names = c("Events", "Participants", "Share of events"),
    omit.coef = "country_iso",
    groups = list(
        "Grievances" = 2:4,
        "Opportunities" = 5:8,
        "Resources" = 9:10
    ),
    caption = "Prais-Winsten regression models of environmental mobilization",
    include.ci = FALSE,
    label = "",
    custom.note = "Robust standard errors in parentheses. %stars",
    custom.gof.rows = list(
        "Country fixed effects" = rep("Yes", 3),
        "Adjusted R sq" = c(summary(M1)$adj.r.squared, summary(M2)$adj.r.squared, summary(M3)$adj.r.squared),
        "Number of obs." = rep(nrow(reg_data), 3)
    ),
    file = paste0(path, "draft/July_2024/Appendix/Tables/reg_model_full_full_sample.tex")
)
# "Lagged dependent variables" = c("YES", "YES", "YES"),
```


```{r, echo = F, message = F, warning = F, fig.height = 8, fig.width = 6, fig.align = "center"}
clean_reg_result <- function(model, DV_name) {
    mod_ci <- cbind(
        coefci(model, level = 0.95, vcov. = vcovPC(model)),
        coefci(model, level = 0.84, vcov. = vcovPC(model)),
        coef(model)
    )

    mod_ci <- as.data.frame(mod_ci)
    mod_ci$term <- rownames(mod_ci)
    colnames(mod_ci)[5] <- "estimate"
    mod_ci$dv_name <- DV_name
    rownames(mod_ci) <- NULL

    return(mod_ci)
}

M1_ci <- clean_reg_result(M1, "Events")
M2_ci <- clean_reg_result(M2, "Participants")
M3_ci <- clean_reg_result(M3, "Share of events")
reg_res <- rbind(M1_ci, M2_ci, M3_ci)

reg_res <- reg_res %>%
    filter(!(grepl("country_iso", term))) %>%
    filter(!(grepl("Intercept", term))) %>%
    filter(!term == "share_demo") %>%
    mutate(term = case_when(
        term == "econ_indic" ~ "Econ. perf.",
        term == "EPI" ~ "Environ. perf.",
        term == "nr_disasters" ~ "Natural disasters",
        term == "gov_env_cmp" ~ "Gov. on env. prot.",
        term == "env_prot_sys" ~ "Parl. on env. prot.",
        term == "green_in_gov" ~ "Green prt. in gov.",
        term == "green_in_parl" ~ "Green prt. in parl.",
        term == "share_organized" ~ "Org. protests",
        term == "NGO_IUCN" ~ "Env. NGOs",
        TRUE ~ term,
    )) %>%
    mutate(term = factor(term, levels = c(
        "Econ. perf.",
        "Environ. perf.", "Natural disasters",
        "Gov. on env. prot.", "Parl. on env. prot.",
        "Green prt. in gov.", "Green prt. in parl.",
        "Org. protests", "Env. NGOs"
    ))) %>%
    mutate(term = fct_rev(term)) %>%
    mutate(type = case_when(
        term %in% c(
            "Econ. perf.", 
            "Environ. perf.", "Natural disasters"
        ) ~ "Grievances",
        term %in% c(
            "Gov. on env. prot.", "Parl. on env. prot.",
            "Green prt. in gov.", "Green prt. in parl."
        ) ~ "Opportunities",
        term %in% c("Org. protests", "Env. NGOs") ~ "Resources"
    )) %>%
    mutate(dv_name = fct_rev(dv_name))

space_between_bars <- 0.65

ggplot(reg_res, aes(y = term)) +
    geom_vline(xintercept = 0, linetype = "dashed", color = "gray60") +
    geom_point(aes(x = estimate, shape = dv_name, color = dv_name), position = position_dodge(width = space_between_bars), size = 2) +
    geom_linerange(aes(
        xmin = `2.5 %`, xmax = `97.5 %`, color = dv_name,
        group = dv_name
    ), size = 0.5, position = position_dodge(width = space_between_bars), key_glyph = "path") +
    geom_linerange(aes(
        xmin = `8 %`, xmax = `92 %`, color = dv_name,
        group = dv_name
    ), size = 1.1, position = position_dodge(width = space_between_bars), key_glyph = "path") +
    facet_wrap(~type, scale = "free_y", ncol = 2) +
    xlab("Prais-Winsten Estimate") +
    scale_color_brewer(name = "", palette = "Set2", guide = guide_legend(reverse = TRUE, nrow = 1, byrow = TRUE)) +
    scale_shape_discrete(name = "", guide = guide_legend(reverse = TRUE, nrow = 2, byrow = TRUE)) +
    theme_linedraw() +
    theme(
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.key.width = unit(1, "cm"),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank()
    )

ggsave(
    plot = last_plot(),
    filename = "coefplot_overtime_full_sample.png",
    path = paste0(path, "draft/July_2024/Appendix/Figures/"),
    width = 7,
    height = 6,
    scale = 0.9,
    dpi = 400
)
```


```{r, echo = F, message = F, warning = F, fig.height = 8, fig.width = 6, fig.align = "center"}

```





